
SUBROUTINE INIT_LAKE(nx,ny,fnmap1,dt)

!The subroutine INIT_LAKE implements
!initialisation of the model (but not initial conditions)

use PHYS_PARAMETERS
use NUMERIC_PARAMS
use PHYS_CONSTANTS
use DRIVING_PARAMS 
use ATMOS
use ARRAYS
use INOUT_PARAMETERS, only : &
& lake_subr_unit_min, &
& lake_subr_unit_max
implicit none

!Input variables
!Reals
real(8), intent(in) :: dt

!Integers
integer(4), intent(in) :: nx
integer(4), intent(in) :: ny

!Characters
character(len=*), intent(in) :: fnmap1

!Output variables

!Local variables
!Reals
real(8), parameter :: hour_sec = 60.*60.
real(8) :: x1

!Integers
integer(4) :: n_unit = lake_subr_unit_min
integer(4) :: i
integer(4) :: j

!Logicals
logical :: uniform_depth

!Characters
character(len=80) :: path2
character(len=80) :: fnmap

data uniform_depth /.true./
data n_unit /999/

!External procedures
character(len=80), external:: fext2

fnmap = fnmap1

call DEFINE_PHYS_PARAMETERS
call DEFINE_PHYS_CONSTANTS
call DEFINE_DRIVING_PARAMS
call ALLOCARRAYS(nx,ny)
call PBLDAT
call COMSOILFORLAKE

if (error_cov==1) then
  if (runmode==2) then
    print*, 'The error covariance estimation algorithm works &
    &only in standalone runs: STOP'
    STOP
  endif
! if (spinup_day/=0) then
!  print*, 'No spinup is allowed when calculating 
!& error covariances: STOP'
!  STOP
! endif 
  if (assim == 1) then
    print*, 'assim = 1 and error_cov = 1: these regimes could not &
    &be turned on simultaneously: STOP'
    STOP
  endif
  if (nx>1.or.ny>1) then
    print*, 'Error covariance calculation is currently adjusted &
    &only for one-point stand-alone runs of the lake model: STOP'
    STOP
  endif
endif
if (assim==1.and.(nx>1.or.ny>1)) then
  print*, 'Data assimilation is currently adjusted &
  &only for one-point stand-alone lake model: STOP'
  STOP
endif

if (dt_out<dt/hour_sec) then
  print*, 'dt_out must be larger or equal to timestep: STOP'
  STOP
endif
   
if (runmode == 2.and.(.not.uniform_depth)) then
  path2 = fext2(fnmap,'dep')
  call CHECK_UNIT(lake_subr_unit_min,lake_subr_unit_max,n_unit)
  open (n_unit, file=path(1:len_trim(path))// &
  & path2(1:len_trim(path2)), status='old')
  call READGRD_LAKE(n_unit,dep_2d,0,nx-1,0,ny-1)
  close (n_unit)
  dep_av = 0.
  x1 = 0.
  do i = 1, nx-1
    do j = 1, ny-1
      if (dep_2d(i,j) > 0.) then
        dep_av = dep_av + dep_2d(i,j)
        x1 = x1 + 1.
      endif
    enddo
  enddo
  dep_av = dep_av / x1
endif

!print*, 'Lake model is initialized'

END SUBROUTINE INIT_LAKE


SUBROUTINE LAKE &
& (tempair1,humair1,pressure1,uwind1,vwind1, &
& longwave1,shortwave1,precip1,Sflux01,zref1,dtl, &
& h10, l10, ls10, hs10, Ts0, Tb0, Tbb0, h_ML0, extwat, extice, &
& kor_in, trib_inflow, Sals0, Salb0, fetch, phi, lam, us0, vs0, &
& Tm, alphax, alphay, a_veg, c_veg, h_veg, area_lake, &
& tsw,hw1,xlew1,cdmw1,ftot, &
& ix,iy,nx,ny,year,month,day,hour,init_T,flag_assim,flag_print, &
& outpath1)

!MAIN SUBROUTINE OF THE MODEL 
!Calculates all the parameters of state of a lake (temperature, currents,
!eddy diffusivity coefficients, snow, ice, soil characteristics, etc.) 
!at the next timestep (i.e. implements transition of variables 
!from i-th time step to (i+1)-th )

!In atmospheric model must be called once each timestep,
!or once each N timesteps, where N = dt_lake/dt_atmos 
!(dt_lake - timestep in the lake model, dt_atmos - timestep the atmospheric model) 

!INPUT VARIABLES:
!tempair1--- air temperature, K;
!humair1 --- specific humidity of air, kg/kg;
!pressure1--- air pressure, Pa;
!uwind1  --- x-component of wind, m/s;
!vwind1  --- y-component of wind, m/s;
!longwave1--- longwave downward radiation, W/m**2;
!shortwave1--- net solar radiation, W/m**2;
!precip1    --- precipitation intensity, m/s;
!xref1   --- the height (m) of level in atmosphere, where temperature, 
!humidity and wind are measured (calculated by atmospheric model);
!trib_inflow1  --- tributary inflow rate, m/s;
!dtl--- timestep (for stability must be not larger than 10-15 min), s;
!ix--- current number of grid point in X-direction;
!iy--- current number of grid point in Y-direction;
!nx--- total number of grid points in X-direction of atmospheric model;
!ny--- total number of grid points in Y-direction of atmospheric model;

!OUTPUT VARIABLES:   
!tsw--- surface temperature of a lake, K;
!hw1--- sensible heat flux from lake, W/m**2; 
!xlew1   --- latent heat flux from lake, W/m**2; 
!cdmw    --- exchange coefficient, m/s ( == wind_speed*exchange_nondimensional_coeficient)

use NUMERIC_PARAMS
use DRIVING_PARAMS
use ATMOS
use PHYS_CONSTANTS
use ARRAYS
!use ASSIM_VAR
use MODI_MASSFLUX_CONVECTION
use CONVECTPAR
use WATER_DENSITY, only: &
& WATER_DENS_TS, &
& DENS_HOSTETLER
use PHYS_FUNC, only: &
& TURB_SCALES, &
& MELTPNT, &
& WATER_FREEZE_MELT, &
& TURB_DENS_FLUX, &
& SINH0, &
& WATER_ALBEDO, &
& EXTINCT_SNOW, &
& UNFRWAT, &
& WI_MAX, &
& WL_MAX
use EVOLUTION_VARIABLES, only : &
& UPDATE_CURRENT_TIMESTEP, &
& UPDATE_NEXT_TIMESTEP
use TURB_CONST, only : &
& niu, lamTM
use METH_OXYG_CONSTANTS, only : &
& ch4_atm0, o2_atm0

implicit none


!------------------------ MAIN VARIABLES ------------------------
!  arrays:    Tw1 and Tw2 - temperature of water, C
! Ti1 and Ti2 - temperature of ice, C
! T - temperature of snow, C
!  functions of time:h1 and h2 - thickness of water, m
! l1 and l2 - thickness of ice, m
! hs1 - thickness of snow, m  
! flag_ice_surf - shows if ice is melting on the top surface, n/d  
!  constants: cw - specific heat of water, J/(kg*K) 
! ci - specific heat of ice, J/(kg*K)
! row0 - density of water, kg/m**3
! roi - density of ice, kg/m**3
! lamw - thermal conductivity of water, J*sec/(m**3*K)
! lami - thermal conductivity of ice, J*sec/(m**3*K)
! L - specific heat of melting, J/kg
! Lwv - specific heat of evaporation, J/kg
! Liv - specific heat of sublimation, J/kg    
! ddz - size of grid element (space), m
! dt - size of grid element (time), sec
! kstratwater - coefficient for linear initial conditions in water
! kstratice - coefficient for linear initial conditions in ice 
!  boundary conditions:
! eFlux - total heat flux on the upper surface,
!   shortwave(1-A)-T**4+longwave-H-LE, J/(m**2*sec)
! Elatent - latent heat flux, J/(m**2*sec) 
! Erad - shortwave radiation, penetrated below a surface, J/(m**2*sec)
! tempair - air temperature (2 m), C
! precip - precipitation, m/sec
!(M) number of layers in water and snow   

!Parameters

integer(4), parameter :: vector_length = 350
real(8), parameter :: year_sec  = 60.*60.*24.*365.
real(8), parameter :: day_sec  = 60.*60.*24.
real(8), parameter :: pi = 3.141592654

!Input variables
!Reals
real(8), intent(in) :: tempair1
real(8), intent(in) :: humair1
real(8), intent(in) :: pressure1
real(8), intent(in) :: uwind1
real(8), intent(in) :: vwind1
real(8), intent(in) :: longwave1
real(8), intent(in) :: zref1
real(8), intent(in) :: shortwave1
real(8), intent(in) :: precip1
real(8), intent(in) :: Sflux01
real(8), intent(in) :: dtl
real(8), intent(in) :: hour

real(8), intent(in) :: h10, l10, ls10, hs10
real(8), intent(in) :: Ts0, Tb0, Tbb0
real(8), intent(in) :: h_ML0
real(8), intent(in) :: extwat, extice
real(8), intent(in) :: kor_in
real(8), intent(in) :: trib_inflow
real(8), intent(in) :: Sals0, Salb0
real(8), intent(in) :: fetch
real(8), intent(in) :: phi, lam
real(8), intent(in) :: us0, vs0
real(8), intent(in) :: Tm
real(8), intent(in) :: alphax, alphay
real(8), intent(in) :: a_veg, c_veg, h_veg
real(8), intent(in) :: area_lake

!Integers
integer(4), intent(in) :: ix, iy
integer(4), intent(in) :: nx, ny
integer(4), intent(in) :: year, month, day
integer(4), intent(in) :: init_T

!Characters
character(len=60), intent(in) :: outpath1

!Logicals
logical, intent(in) :: flag_assim
logical, intent(in) :: flag_print

!Output variables
!Reals
real(8), intent(out) :: hw1, xlew1, cdmw1, tsw
real(8), intent(out) :: ftot

!Local variables
!Reals
real(8) :: l2, h2, ls2
real(8) :: dhwhigh, dhwlow, dhwp, dhwe
real(8) :: dhihigh, dhilow, dhip, dhis, dhif
real(8) :: totalevap, totalprecip, totalmelt, totalpen, totalwat
real(8) :: snowmass, snowmass_init
real(8) :: totalerad, totalhflux
real(8) :: Tf_old1
real(8) :: h_ML0zv
real(8) :: dhwls
real(8) :: x, xx
real(8) :: kor
real(8) :: Ti10 ! Initial temperature of the ice surface (if l10 > 0)

real(8), dimension(1:vector_length) :: a, b, c, d 

real(8) :: lams(ms) 
real(8) :: dz
real(8) :: ST
real(8) :: PGR
real(8) :: TGROLD
real(8) :: QGROLD
real(8) :: RADIAT
real(8) :: WSOLD
real(8) :: SNOLD
real(8) :: ZS
real(8) :: thsoil
real(8) :: whsoil
real(8) :: bOLD
real(8) :: RF1, RF2
real(8) :: SNMELT
real(8) :: HSNOW
real(8) :: HS
real(8) :: ES
real(8) :: TGRNEW
real(8) :: QGRNEW
real(8) :: WSNEW
real(8) :: SNNEW
real(8) :: RUNOFF
real(8) :: ElatOld
real(8) :: HFold
real(8) :: PRSold
real(8) :: extinct(ms)

real(8) :: AL
real(8) :: DLT
real(8) :: DVT
real(8) :: ALLL
real(8) :: DL
real(8) :: ALV
real(8) :: DV
real(8) :: Z
real(8) :: T
real(8) :: WL
real(8) :: WV
real(8) :: WI
real(8) :: dens

real(8) :: dt
real(8) :: q(110)
real(8) :: roughness
real(8) :: emissivity
real(8) :: albedo
real(8) :: aM, bM
real(8) :: relhums
real(8) :: b0
real(8) :: tau_air, tau_i, tau_gr
real(8) :: eflux0_kinem_water

real(8), dimension(1:vector_length) :: Temp

real(8) :: flux1, flux2

real(8), external :: DZETA
real(8), external :: DZETAI


!Integers    
integer(4) :: i, j
integer(4) :: i_ML
integer(4) :: flag_snow
integer(4) :: itop
integer(4) :: nstep_meas
integer(4) :: n_1cm, n_5cm
integer(4) :: layer_case
integer(4) :: ndec

!Logicals
logical :: soil_ts_ext_iter
logical :: uniform_depth
logical :: flag
logical :: accum = .false.
logical :: add_to_winter

!Characters
character(len=60) :: outpath

common /watericesnowarr/ lams,q
common /SOILDAT/ dz(ms),itop
common /bL/ ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD, &
& ZS,thsoil,whsoil,bOLD,RF1,RF2,SNMELT,HSNOW, &
& HS,ES,TGRNEW,QGRNEW,WSNEW,SNNEW,RUNOFF, &
& ElatOld,HFold,PRSold,extinct
common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
& ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML), &
& dens(ms)
common /ts_ext_iter/ Tf_old1,soil_ts_ext_iter 
common /surface/ roughness,emissivity,albedo,aM,bM,relhums
common /out/ outpath

data uniform_depth /.true./
data nstep_meas /0/
 
SAVE
  
!BODY OF PROGRAM
    
tempair = tempair1-273.15
humair = humair1
pressure = pressure1
!pressure = 1.d+5
uwind = uwind1
vwind = vwind1
zref  = zref1
precip = precip1
Sflux0 = Sflux01*(roa0/row0)
outpath = outpath1

if (kor_in == -999.d0) then
  kor = 2.d0*omega*dsin(phi*pi/180.d0)
else
  kor = kor_in
endif

if(ifrad == 1) then
  longwave  = longwave1
  shortwave = shortwave1
elseif (ifrad == 0) then
  longwave  = 0.d0
  shortwave = 0.d0
endif

if (uwind == 0) uwind = 0.1
if (vwind == 0) vwind = 0.1
wind = dsqrt(uwind**2+vwind**2)
!wind10 = wind*dlog(10./roughness0)/dlog(zref/roughness0)

!Control of the input atmospheric forcing
if (tempair>60.or.tempair<-90) then
  print*, 'The air temperature ', tempair, 'is illegal: STOP'
  STOP
elseif (dabs(humair)>1.) then
  print*, 'The air humidity ', humair, 'is illegal: STOP'
  STOP
elseif (pressure > 110000 .or. pressure < 90000.) then
  print*, 'The air pressure ', pressure, 'is illegal: STOP'
  STOP
elseif (dabs(uwind)>200.) then
  print*, 'The x-component of wind ', uwind, 'is illegal: STOP'
  STOP
elseif (dabs(vwind)>200.) then
  print*, 'The y-component of wind ', vwind, 'is illegal: STOP'
  STOP
elseif (dabs(longwave)>1000.) then
  print*, 'The longwave radiation ',longwave,'is illegal: STOP'
  STOP
elseif (dabs(shortwave)>1400.) then
  print*, 'The shortwave radiation ', shortwave, 'is illegal: STOP'
  STOP
elseif (dabs(precip)>1.) then
  print*, 'The atmospheric precipitation ',precip, &
  & 'is illegal: STOP'
  STOP
endif

dt = dtl
dt05 = 0.5d0*dt
dt_keps = dt/real(nstep_keps)
dt_inv = 1.d0/dt
dt_inv05 = 0.5d0*dt_inv

if (init(ix,iy)==0) then
! Specification of initial profiles   

!  write(*,'(a)') 'The thickness of (1) dzeta and (2) z layers in water coloumn'
!  write(*,'(a,i5,a,i5)') 'in the point ix = ', ix, ', iy = ', iy
!  do i = 1, M
!    write(*,'(a5,i4,a7,f10.5,a7,f10.5)') ' N ', i, ' (1) ', ddz(i), ' (2) ', ddz(i)*h10
!  enddo
  
  rosoil(:) = 1200. ! Bulk soil density for initialization
  call INIT_VAR &
  ( M, Mice, ns, ms, ml, Tinitlength, init_T, skin, zero_model, &
  & h10, l10, hs10, ls10, tempair, Ts0, Tb0, Tm, Tbb0, &
  & Sals0, Salb0, us0, vs0, h_ML0, &
  & rosoil, por, &
  & zTinitprof, Tinitprof, &
  
  & flag_snow, flag_snow_init, itop, nstep, &
  & h1, l1, hs1, ls1, &
  & veg, snmelt, cdmw2, velfrict_prev, &
  & roughness, eflux0_kinem, Elatent, &
  & totalevap, totalmelt, totalprecip, totalwat, totalpen, &
  & time, dhwfsoil, dhw, dhw0, dhi, dhi0, dls0, &
  
  & E1, eps1, Tsoil1, wi1, wl1, Sals1, &
  & rootss, qsoil, TgrAnn, qwater(1,1), &
  & oxyg(1,1), carbdi(1,1), Sal1, u1, v1, Tw1, Tskin, T_0dim, &
  & Ti1, Tis1, z_full, ddz, &
  & dz, T, wl, dens, lamw, &
  & dzeta_int, zsoil, pressure )    
      
  if (runmode==2.and.(uniform_depth .eqv. .false.)) then
    if (ix<=nx-1.and.iy<=ny-1) then
      h1 = dep_2d(ix,iy)
    else
      h1 = dep_av
    endif  
    if (h1<=0) then
      print*, 'Negative or zero initial depth at ix=',ix, &
      & 'iy=',iy,'h1=',h1,':terminating program'
      STOP 
    endif
  endif

! if (assim == 1.or.error_cov == 1) then
!   call measured_data_reader
! endif

else

call UPDATE_CURRENT_TIMESTEP ( &
& ix, iy, nx, ny, M, Mice, ns, ms, ml, &
& l1, h1, ls1, hs1, &
& u1, v1, &
& E1, eps1, k_turb_T_flux, &
& Tsoil1, Sals1, wi1, wl1, &
& Tw1, Sal1, lamw, &
& Tskin(1), &
& Ti1, Tis1, &
& dz, T, wl, dens, &
& qwater(1,1), qsoil, &
& oxyg(1,1), carbdi(1,1), &
& snmelt, &
& cdmw2, &
& time, &
& dhwfsoil, &
& Elatent, &
& dhw, &
& dhw0, &
& dhi, &
& dhi0, &
& dls0, &
& velfrict_prev, &
& roughness, &
& eflux0_kinem, &
& tot_ice_meth_bubbles, &

& flag_snow, flag_snow_init, &
& itop, &
& nstep )

endif

! Linear decrease of horizontal cross-section with depth
!do i = 1, M+1
!  area_int(i) = (1.-DZETA(real(i))*0.9 )*area_lake
!enddo
!do i = 1, M
!  area_half(i) = (1.-DZETA(real(i)+0.5)*0.9 )*area_lake
!enddo

! Constant cross-section 
do i = 1, M+1
  area_int(i) = area_lake
enddo
do i = 1, M
  area_half(i) = area_lake
enddo

!if (flag_assim==.true.) then
!Correcting model fields at the current timestep 
!using least-square data assimilation   
! do i = 1, M+1
!  j = (i-1)*n_modvar
!  x_mod(1,j+1)= Tw1(i)
!  x_mod(1,j+2)= E1(i)
! enddo
! call assim_lst_sq 
!endif
   
time = time + dt
nstep = nstep + 1

layer_case = 1
if (h1 > 0  .and. l1 > 0) layer_case = 2 
if (h1 ==0  .and. l1 > 0) layer_case = 3
if (h1 ==0  .and. l1 ==0) layer_case = 4

if (layer_case == 1 .and. &
& WATER_FREEZE_MELT(Tw1(1), 0.5*ddz(1)*h1, Meltpnt(Sal1(1)), +1) .and. &
& h1 - min_ice_thick * roi/row0 > min_water_thick) then
  layer_case = 2
endif 

if (layer_case == 3 .and. &
& WATER_FREEZE_MELT(Ti1(Mice+1), 0.5*ddzi(Mice)*l1, Meltpnt(0.d0), -1) .and. &
& l1 - min_water_thick * row0_d_roi > min_ice_thick) then
  layer_case = 2
endif  

if (layer_case == 3 .and. &
& WATER_FREEZE_MELT(Ti1(1), 0.5*ddzi(1)*l1, Meltpnt(0.d0), -1) .and. &
& l1 - min_water_thick * row0_d_roi > min_ice_thick .and. flag_snow == 0) then
  h1 = min_water_thick
  Tw1 = Meltpnt(0.d0) + T_phase_threshold 
  Sal1 = 1.d-5
  ls1 = l1 - min_water_thick * row0_d_roi
  Tis1 = Ti1
  l2 = 0
  Ti2 = 0
  if (ls1<0.009999) print*, 'Thin layer! - ls'
  layer_case = 1 
endif

  
if1: IF (layer_case == 1) THEN 

!---------------------------- CASE 1: LIQUID WATER LAYER AT THE TOP ("Summer")-----------------------

!Check the bottom layer water temperature
if (WATER_FREEZE_MELT(Tw1(M+1), 0.5*ddz(M)*h1, Meltpnt(Sal1(M+1)), +1) .and. &
& ls1 == 0 .and. h1 - min_ice_thick * roi/row0 > min_water_thick) then
  ls1 = min_ice_thick
  Tis1 = Meltpnt(Sal1(M+1)) - T_phase_threshold
  h1 = h1 - min_ice_thick*roi/row0
endif   

!Check the liquid temperature to be above melting point
do i = 2, M
  if (Tw1(i)<Meltpnt(Sal1(i))) Tw1(i)=Meltpnt(Sal1(i))
enddo

!Calculation of water density profile at the current timestep,
!taking into account only temperature
do i = 1, M+1
  row(i) = WATER_DENS_TS(Tw1(i),Sal1(i)) !0.d0
!  row(i) = DENS_HOSTETLER(Tw1(i))
enddo

!The salinity and salinity flux at the surface are passed to
!TURB_DENS_FLUX as zeros since they are currently
!                                                      ______
!not taken into account in calculation of density flux w'row'	
 
turb_density_flux = TURB_DENS_FLUX(eflux0_kinem,0.d0,Tw1(1),0.d0)
Buoyancy0 = g/row0*turb_density_flux

do i = 1, M+1
  z_full  (i) = dzeta_int(i)*h1
enddo

if (massflux == 1) then
!Updating density profile for massflux convection
!Updating z-grid properties 
  do i = 1, M
    dz_full (i) = ddz(i)  *h1
    z_half  (i) = dzeta_05int(i)*h1
  enddo
  call MASSFLUX_CONVECTION( &
  & nstep,dt,z_half,z_full,dz_full, &
  & turb_density_flux,eflux0_kinem, &
! The latent heat flux and salinity flux at the surface are passed as zeros,
! as far as currently they are not taken 
! into account in massflux conection parameterization
  & 0.d0,0.d0, &
! The surface wind stress components are passed as zeros,
! as far as currently they are not taken 
! into account in massflux conection parameterization
  & 0.d0,0.d0, &
  & u1,v1,w1,E1(1:M), &
  & u1,v1,w1,row,Tw1,Sal1, &
  & PEMF,PDENS_DOWN,PT_DOWN,PSAL_DOWN,zinv,1,M)
! Debugging purposes only
! PEMF = 0.d0
! call TEMP_MASSFLUX(PEMF,PT_DOWN,ddz,h1,dt,Tw2,T_massflux)
  do i = 1, M
    pt_down_f(i) = 0.5d0*(PT_DOWN(i) + PT_DOWN(i+1))
  enddo
endif

! The solar radiation, absorbed by the water surface
if (varalb == 0) then
  Erad = shortwave*(1-albedoofwater)
elseif (varalb == 1) then
  Erad = shortwave* &
  & (1-WATER_ALBEDO( SINH0(year,month,day,hour,phi) ) ) ! Note: the time here must be a local one
endif  ! time, not UTC

!Radiation fluxes in layers
do i = 1, M
  SR(i)=Erad*(1-sabs)*dexp(-extwat*dzeta_05int(i)*h1)
enddo
if (ls1 /= 0.d0) then
  do i = 1, Mice
    SRdi(i) = Erad*(1-sabs)*dexp(-extwat*h1)* &
    & (1-albedoofice)*dexp(-extice*dzetai_05int(i)*ls1)
  enddo
endif

!Calculation of the layer's thickness increments
if (ls1==0) then
  ice=0;snow=0;water=1;deepice=0
  dhwp = precip*dt
  dhwe = - Elatent/(Lwv*row0)*dt
  dhw = dhwp + dhwe + dhwfsoil
  dhw0 = dhwp + dhwe
  dls = 0.
else
  ice=0;snow=0;water=1;deepice=1
  flux1 = -lamw(M)*(Tw1(M+1)-Tw1(M))/(ddz(M)*h1) + SR(M) + &
  & cw_m_row0*(dhw-dhw0)*(Tw1(M+1)-Tw1(M))/(2.*dt) + &
  & cw_m_row0*PEMF(M)*(pt_down_f(M)-0.5d0*(Tw1(M+1)+Tw1(M)) )
  flux2 = -lami*(Tis1(2)-Tis1(1))/(ddzi(1)*ls1) + SRdi(1) + &
  & ci_m_roi*dls0*(Tis1(2)-Tis1(1))/(2.*dt)
  dhwls = dt*(flux1 - flux2)/(row0_m_Lwi)
!  dhwls = -lamw(M)*dt/(ddz(M)*h1*row0_m_Lwi)*(Tw1(M+1) - Tw1(M))+
!&  lami*dt/(ddz(1)*ls1*row0_m_Lwi)*(Tis1(2) - Tis1(1))+
!&  (1-albedoofice)*Erad*(1-sabs)*dexp(-extwat*h1)/(row0_m_Lwi)*dt 
  dls = -dhwls*row0_d_roi
  dls0 = dls
  dhwp = precip*dt
  dhwe = - Elatent/(Lwv*row0)*dt
  dhw = dhwp + dhwe + dhwfsoil + dhwls
  dhw0 = dhwp + dhwe
endif

!Calculation of water current's velocities and turbulent characteristics
if (Turbpar /= 8) then
  call MOMENT_SOLVER(ix, iy, nx, ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, &
  & alphax, alphay, dt, b0, tau_air, tau_i, tau_gr)
endif

if (Turbpar == 2) then
  do i = 1, nstep_keps
    call K_EPSILON(ix,iy,nx,ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, dt_keps, &
    & b0, tau_air, tau_i, tau_gr, roughness)
  enddo
else
! Calculation of the of eddy difusivity for temperature
  SELECT CASE (Turbpar)
    CASE (8)
!     Hostetler formulation
      call ED_TEMP_HOSTETLER &
      & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
      KT(:) = cw_m_row0*KT(:)
    CASE (1)
      do i = 1, M
        KT(i) = (k2(i)-niu)*cw_m_row0 !*0.01
      enddo
    CASE DEFAULT
      KT = lamTM*k2*cw_m_row0 !*0.01
  ENDSELECT 
endif
  
do i = 1, M
  lamw(i) = lamw0 + KT(i) !+KC(i) !+KLengT(i)
enddo
lamsal(:) = lamw(:)/cw_m_row0 * alsal
lammeth(:) = lamw(:)/cw_m_row0 * almeth
lamoxyg(:) = lamw(:)/cw_m_row0 * aloxyg
lamcarbdi(:) = lamw(:)/cw_m_row0 * alcarbdi

!Calculation of the whole temperature profile
call SOIL_COND_HEAT_COEF()
CALL S_DIFF(dt)
CALL T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
& extwat, extice, fetch, dt)
call SOILFORLAKE &
& (ix,iy,nx,ny,year,month,day,hour,phi, &
& extwat,extice,fetch,dt,a,b,c,d,Temp)
call METHANE &
& (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil,TgrAnn, &
& ddz, Tw2, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
& fplant, febul, fdiff, ftot, fdiff_lake_surf, &
& plant_sum,bull_sum,oxid_sum,rprod_sum, &
& anox,M,ns,dt,rprod_total_oldC,rprod_total_newC,h_talik,tot_ice_meth_bubbles,add_to_winter)
qmethane(1:M+1) = qwater(1:M+1,2)
qmethane(M+2:M+ns) = qsoil(2:ns)
!call METHANE2 & ! two-meth
!& (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!& ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
!& fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!& plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!& anox,M,ns,dt,rprod_total_oldC,rprod_total_newC, & ! two-meth
!& h_talik,tot_ice_meth_bubbles2) ! two-meth
!qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
call OXYGEN (M, Tw2, lamoxyg, ddz, pressure, wind10, h1, l1, ls1, dhw, dhw0, dt, oxyg)
call CARBON_DIOXIDE (M, Tw2, lamcarbdi, ddz, pressure, wind10, &
& h1, l1, ls1, dhw, dhw0, dt, carbdi)
call METHANE_OXIDATION(M, dt, oxyg(1,2), qwater(1,2), carbdi(1,2), Tw2)
!call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), carbdi, Tw2) ! two-meth


!Diagnostic calculations
do i = 1, M
  k_turb_T_flux(i) = - lamw(i)/(cw_m_row0)*( Tw2(i+1) - Tw2(i) )/ &
  & (ddz(i)*h1)
  T_massflux(i) = PEMF(i)*(pt_down_f(i)-0.5d0*(Tw2(i+1)+Tw2(i)))
enddo

!Adding the heat input from tributaries
if (tribheat==1) call TRIBTEMP(dt,area_lake,Tw2)

!The scales of turbulence
call TURB_SCALES(k_turb_T_flux, T_massflux, eflux0_kinem, h1, &
& turb_density_flux, Buoyancy0, H_mixed_layer, w_conv_scale, &
& T_conv_scale)

if (zero_model == 1) then
! Note that zero model is impemented currently only for
! 1. one-point simulation
! 2. open water conditions
  call ZERODIM_MODEL &
  & (h1, dt, &
  & shortwave, longwave, tempair, humair, pressure, &
  & uwind, vwind, zref, T_0dim)
endif
  
h2 = h1 + dhw
ls2 = ls1 + dls
l2 = 0.

ENDIF if1   

if2: IF (layer_case==2) THEN

!---------------------------CASE 2: WATER AND ICE ("Winter")------------------------------


! Creation of the initial thin ice layer
if (l1 == 0.) then
  l1 = min_ice_thick
  h1 = h1 - min_ice_thick * roi/row0
  if (tempair<0.) then
    do i=1, Mice+1
      Ti1(i) = tempair + float(i-1)/float(Mice) * &
      & (Meltpnt(Sal1(1))-tempair)
    enddo 
  else 
    Ti1 = Meltpnt(Sal1(1)) - T_phase_threshold
  endif
end if 

!Check for water temperature to be higher than melting point
do i=2,M
  if (Tw1(i)<Meltpnt(Sal1(i))) Tw1(i)=Meltpnt(Sal1(i))
enddo

!Creation of the initial thin water layer
if (h1 == 0.) then
  h1 = min_water_thick
  Sal1 = 1.d-5
  Tw1 = Meltpnt(0.d0) + T_phase_threshold
  l1 = l1 - min_water_thick * row0_d_roi
end if

!Creation of the initial thin bottom ice layer
if (WATER_FREEZE_MELT(Tw1(M+1), 0.5*ddz(M)*h1, Meltpnt(Sal1(M+1)), +1) .and. &
& ls1 == 0 .and. h1 - min_ice_thick * roi/row0 > min_water_thick) then
  ls1 = min_ice_thick
  Tis1 = Meltpnt(Sal1(M+1)) - T_phase_threshold
  h1 = h1 - min_ice_thick * roi/row0
endif

!Calculation of water density profile at the current timestep,
!taking into account only temperature
do i = 1, M+1
  row(i) = WATER_DENS_TS(Tw1(i),Sal1(i)) !0.d0
!  row(i) = DENS_HOSTETLER(Tw1(i))
enddo

! Note: the salinity effects and partial ice coverage 
! are currently neglected in density flux at the ice-water interface
eflux0_kinem_water = - lamw(1)*(Tw1(2)-Tw1(1))/(ddz(1)*h1*cw_m_row0) ! The first-order approximation
turb_density_flux = TURB_DENS_FLUX(eflux0_kinem_water,0.d0,Tw1(1),0.d0)
Buoyancy0 = g/row0*turb_density_flux

do i = 1, M+1
  z_full  (i) = dzeta_int(i)*h1
enddo

!Currently the EDMF parameterization of convection is not implemented in the code
!during 'winter' (if the ice layer is present over
!water layer). This seems to be reasonable, since the temperature profile
!is stable in this case - until spring convection under thin ice
!develops.

PEMF = 0.d0
pt_down_f = 0.d0
  
!CASE 2.1: WATER, ICE AND SNOW 

if (flag_snow == 1) then

!Radiation fluxes in layers
  SR_botsnow = shortwave*(1-albedoofsnow)
  if (flag_snow_init /= 1) then
    do i = itop, ms-1
      SR_botsnow = SR_botsnow*EXTINCT_SNOW(dens(i))**dz(i)
    enddo
  endif
  do i = 1, Mice
    SRi(i) = SR_botsnow*dexp(-extice*dzetai_05int(i)*l1)
  enddo
  do i = 1, M
    SR(i) = SR_botsnow*dexp(-extice*l1)* &
    & dexp(-extwat*dzeta_05int(i)*h1)
  enddo
  if (ls1 /= 0.d0) then
    do i = 1, Mice
      SRdi(i) = SR_botsnow*dexp(-extice*l1)*dexp(-extwat*h1)* &
      & (1-albedoofice)*dexp(-extice*dzetai_05int(i)*ls1)
    enddo
  endif


! Calculation of the layer's thickness increments
  flux1=-lami*(Ti1(Mice+1)-Ti1(Mice))/(ddzi(Mice)*l1)+SRi(Mice) + &
  & ci_m_roi*(dhi-dhi0)*(Ti1(Mice+1)-Ti1(Mice))/(2.d0*dt)
  flux2 = -lamw(1)*(Tw1(2)-Tw1(1))/(ddz(1)*h1) + SR(1) + &
  & cw_m_row0*dhw0*(Tw1(2)-Tw1(1))/(2.*dt) + &
  & cw_m_row0*PEMF(1)*(pt_down_f(1)-0.5d0*(Tw1(2)+Tw1(1)) )
  dhwlow = dt*(flux1 - flux2)/(row0_m_Lwi)
  dhilow = -dhwlow*row0_d_roi
  if (Ti1(1) > Meltpnt(0.d0) + T_phase_threshold) then
    dhwhigh = (Ti1(1) - Meltpnt(0.d0) - T_phase_threshold) * &
    & ci_m_roi*l1*ddzi(1)*0.5d0/(row0_m_Lwi)
    dhihigh = -dhwhigh*row0_d_roi
  else
    dhwhigh = 0.d0
    dhihigh = 0.d0
  endif
  if (ls1/=0) then
    ice=1;snow=1;water=1;deepice=1
    flux1 = -lamw(M)*(Tw1(M+1)-Tw1(M))/(ddz(M)*h1) + SR(M) + &
    & cw_m_row0*(dhw-dhw0)*(Tw1(M+1)-Tw1(M))/(2.d0*dt) + &
    & cw_m_row0*PEMF(M)*(pt_down_f(M)-0.5d0*(Tw1(M+1)+Tw1(M)) )
    flux2 = -lami*(Tis1(2)-Tis1(1))/(ddzi(1)*ls1) + SRdi(1) + &
    & ci_m_roi*dls0*(Tis1(2)-Tis1(1))/(2.d0*dt)
    dhwls = dt*(flux1 - flux2)/(row0_m_Lwi)
    dls = -dhwls*row0_d_roi
    dhw = dhwhigh + dhwlow + snmelt*dt + dhwfsoil + dhwls
    dhw0 = dhwlow + dhwhigh + snmelt*dt
    dls0 = dls
  else
    ice=1;snow=1;water=1;deepice=0
    dls = 0.
    dhwls = 0.
    dhw = dhwhigh + dhwlow + snmelt*dt + dhwfsoil
    dhw0 = dhwlow + dhwhigh + snmelt*dt
  endif
  dhi = dhilow + dhihigh 
  dhi0 = dhihigh

!Calculation of water current's velocities and turbulent characteristics
  if (Turbpar /= 1) then
  
    if (Turbpar /= 8) then
      call MOMENT_SOLVER(ix, iy, nx, ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, &
      & alphax, alphay, dt, b0, tau_air, tau_i, tau_gr)  
    endif
    
    if (Turbpar == 2) then
      do i = 1, nstep_keps
        call K_EPSILON(ix,iy,nx,ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, dt_keps, &
        & b0, tau_air, tau_i, tau_gr, roughness)
      enddo
    else
!     Calculation of the of eddy difusivity for temperature
      SELECT CASE (Turbpar)
        CASE (8)
!         Hostetler formulation
!          call ED_TEMP_HOSTETLER &
!          & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
          KT(:) = 0.
        CASE (1)
          do i = 1, M
            KT(i) = (k2(i)-niu)*cw_m_row0 !*0.01
          enddo
        CASE DEFAULT
          KT = lamTM*k2*cw_m_row0 !*0.01
      ENDSELECT 
    endif
    
    do i = 1, M
      lamw(i) = lamw0 + KT(i) !+KC(i)
    enddo  
    
  else
    lamw = lamw0*3.d0
  endif
  
  lamsal(:) = lamw(:)/cw_m_row0 * alsal
  lammeth(:) = lamw(:)/cw_m_row0 * almeth
  lamoxyg(:)= lamw(:)/cw_m_row0 * aloxyg
  lamcarbdi(:) = lamw(:)/cw_m_row0 * alcarbdi

!Calculation of the whole temperature profile 
  call SNOW_COND_HEAT_COEF()
  call SOIL_COND_HEAT_COEF()
  call S_DIFF(dt)
  call T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat, extice, fetch, dt)
  call SOILFORLAKE( &
  & ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat,extice,fetch,dt,a,b,c,d,Temp)
  call SNOWTEMP(ix,iy,nx,ny,year,month,day,hour,snowmass, &
  & snowmass_init,a,b,c,d,Temp,phi,extwat,extice,fetch,dt)
  call METHANE &
  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil,TgrAnn, &
  & ddz, Tw2, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
  & fplant, febul, fdiff, ftot, fdiff_lake_surf, &
  & plant_sum,bull_sum,oxid_sum,rprod_sum, &
  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC,h_talik,tot_ice_meth_bubbles,add_to_winter)
  qmethane(1:M+1) = qwater(1:M+1,2)
  qmethane(M+2:M+ns) = qsoil(2:ns)
!  call METHANE2 & ! two-meth
!  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
!  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC, & ! two-meth
!  & h_talik,tot_ice_meth_bubbles2) ! two-meth
!  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
  call OXYGEN (M, Tw2, lamoxyg, ddz, pressure, wind10, h1, l1, ls1, dhw, dhw0, dt, oxyg)
  call CARBON_DIOXIDE (M, Tw2, lamcarbdi, ddz, pressure, wind10, &
  & h1, l1, ls1, dhw, dhw0, dt, carbdi)
  call METHANE_OXIDATION(M, dt, oxyg(1,2), qwater(1,2), carbdi(1,2), Tw2)
!  call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), carbdi, Tw2) ! two-meth
  
  h2 = h1 + dhw
  l2 = l1 + dhi 
  ls2 = ls1 + dls
  
else

!CASE 2.2: WATER AND ICE WITHOUT SNOW

!Radiation fluxes in layers
  Erad = shortwave*(1-albedoofice)
  do i = 1, Mice
    SRi(i) = Erad*dexp(-extice*dzetai_05int(i)*l1)
  enddo
  do i = 1, M
    SR(i) = Erad*dexp(-extice*l1)* &
    & dexp(-extwat*dzeta_05int(i)*h1)
  enddo
  if (ls1 /= 0.d0) then
    do i = 1, Mice
      SRdi(i) = Erad*dexp(-extice*l1)*dexp(-extwat*h1)* &
      & (1-albedoofice)*dexp(-extice*dzetai_05int(i)*ls1)
   enddo
  endif

!Calculation of the layer's thickness increments
  if (Ti1(1) > Meltpnt(0.d0) + T_phase_threshold) then
    dhwhigh = (Ti1(1) - Meltpnt(0.d0) - T_phase_threshold) * &
    & ci_m_roi*ddzi(1)*l1*0.5d0/(row0_m_Lwi)
    dhis = 0. 
    dhihigh = -dhwhigh*row0_d_roi
  endif
! Creation of the initial thin snow layer - to be considered at the next timestep
  if (precip>0.and.tempair<0.and.l1>0.05) then
    flag_snow = 1
    hs1 = 0.02
    itop = ms-2
    dhip = 0
  else
    dhip = precip*row0_d_roi*dt
  end if
  flux1=-lami*(Ti1(Mice+1)-Ti1(Mice))/(ddzi(Mice)*l1)+SRi(Mice) + &
  &  ci_m_roi*(dhi-dhi0)*(Ti1(Mice+1)-Ti1(Mice))/(2.d0*dt)
  flux2 = -lamw(1)*(Tw1(2)-Tw1(1))/(ddz(1)*h1) + SR(1) + &
  &  cw_m_row0*dhw0*(Tw1(2)-Tw1(1))/(2.*dt) + &
  &  cw_m_row0*PEMF(1)*(pt_down_f(1)-0.5d0*(Tw1(2)+Tw1(1)) )
  dhwlow = dt*(flux1 - flux2)/(row0_m_Lwi)
  dhilow = -dhwlow*row0_d_roi

  if (Ti1(1) < Meltpnt(0.d0) + T_phase_threshold) then
    dhis = - Elatent/(roi*Liv)*dt
    dhwhigh = 0
    dhihigh = 0
  endif
  dhi = dhihigh + dhilow + dhis + dhip
  dhi0 = dhis + dhihigh + dhip
  if (ls1/=0) then
    ice=1;snow=0;water=1;deepice=1 
    flux1 = -lamw(M)*(Tw1(M+1)-Tw1(M))/(ddz(M)*h1) + SR(M) + &
    & cw_m_row0*(dhw-dhw0)*(Tw1(M+1)-Tw1(M))/(2.d0*dt) + &
    & cw_m_row0*PEMF(M)*(pt_down_f(M)-0.5d0*(Tw1(M+1)+Tw1(M)) )
    flux2 = -lami*(Tis1(2)-Tis1(1))/(ddzi(1)*ls1) + SRdi(1) + &
    & ci_m_roi*dls0*(Tis1(2)-Tis1(1))/(2.d0*dt)
    dhwls = dt*(flux1 - flux2)/(row0_m_Lwi)
    dls = -dhwls*row0_d_roi
    dls0 = dls
    dhw = dhwhigh + dhwfsoil + dhwls + dhwlow
    dhw0 = dhwhigh+dhwlow
  else
    ice=1;snow=0;water=1;deepice=0
    dhw = dhwhigh + dhwfsoil + dhwlow
    dhw0 = dhwhigh + dhwlow
    dls = 0.
  endif

!Calculation of water current's velocities and turbulent characteristics
  if (Turbpar/=1) then
  
    if (Turbpar /= 8) then
      call MOMENT_SOLVER(ix, iy, nx, ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, &
      & alphax, alphay, dt, b0, tau_air, tau_i, tau_gr)  
    endif
    
    if (Turbpar == 2) then
      do i = 1, nstep_keps
        call K_EPSILON(ix,iy,nx,ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, dt_keps, &
        & b0, tau_air, tau_i, tau_gr, roughness)
      enddo
    else
!     Calculation of the of eddy difusivity for temperature
      SELECT CASE (Turbpar)
        CASE (8)
!       Hostetler formulation
!          call ED_TEMP_HOSTETLER &
!          & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
          KT(:) = 0.
        CASE (1)
          do i = 1, M
            KT(i) = (k2(i)-niu)*cw_m_row0 !*0.01
          enddo
        CASE DEFAULT
          KT = lamTM*k2*cw_m_row0 !*0.01
      ENDSELECT 
    endif
    
    do i = 1, M
      lamw(i) = lamw0 + KT(i) !+KC(i)
    enddo 
    
  else
    lamw = lamw0*3.d0
  endif
  
  lamsal(:) = lamw(:)/cw_m_row0 * alsal
  lammeth(:) = lamw(:)/cw_m_row0 * almeth
  lamoxyg(:) = lamw(:)/cw_m_row0 * aloxyg
  lamcarbdi(:) = lamw(:)/cw_m_row0 * alcarbdi

! Calculation of the whole temperature profile
  call SOIL_COND_HEAT_COEF()
  call S_DIFF(dt)
  call T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat, extice, fetch, dt)
  call SOILFORLAKE( &
  & ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat,extice,fetch,dt,a,b,c,d,Temp)
  call METHANE &
  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil,TgrAnn, &
  & ddz, Tw2, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
  & fplant, febul, fdiff, ftot, fdiff_lake_surf, &
  & plant_sum,bull_sum,oxid_sum,rprod_sum, &
  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC,h_talik,tot_ice_meth_bubbles,add_to_winter)
  qmethane(1:M+1) = qwater(1:M+1,2)
  qmethane(M+2:M+ns) = qsoil(2:ns)
!  call METHANE2 & ! two-meth
!  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
!  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC, & ! two-meth
!  & h_talik,tot_ice_meth_bubbles2) ! two-meth
!  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
  call OXYGEN (M, Tw2, lamoxyg, ddz, pressure, wind10, h1, l1, ls1, dhw, dhw0, dt, oxyg)
  call CARBON_DIOXIDE (M, Tw2, lamcarbdi, ddz, pressure, wind10, &
  & h1, l1, ls1, dhw, dhw0, dt, carbdi)
  call METHANE_OXIDATION(M, dt, oxyg(1,2), qwater(1,2), carbdi(1,2), Tw2)
!  call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), carbdi, Tw2) ! two-meth

  h2 = h1 + dhw
  l2 = l1 + dhi
  ls2 = ls1 + dls
 
endif

!Diagnostic calculations
do i = 1, M
  k_turb_T_flux(i) = - lamw(i)/(cw_m_row0)*( Tw2(i+1) - Tw2(i) )/ &
  & (ddz(i)*h1)
  T_massflux(i) = PEMF(i)*(pt_down_f(i)-0.5d0*(Tw2(i+1)+Tw2(i)))
enddo

if (tribheat==1) call TRIBTEMP(dt,area_lake,Tw2)

!The scales of turbulence
call TURB_SCALES(k_turb_T_flux, T_massflux, eflux0_kinem_water, h1, &
& turb_density_flux, Buoyancy0, H_mixed_layer, w_conv_scale, &
& T_conv_scale)
 
ENDIF if2

!CASE 3: ICE WITHOUT WATER

if3: IF (layer_case==3) THEN
    
!CASE 3.1: ICE WITH SNOW

if (flag_snow == 1) then

! Radiation fluxes
  SR_botsnow = shortwave*(1-albedoofsnow)
  if (flag_snow_init /= 1) then
    do i = itop, ms-1
      SR_botsnow = SR_botsnow*EXTINCT_SNOW(dens(i))**dz(i)
    enddo
  endif
  do i = 1, Mice
    SRi(i) = SR_botsnow*dexp(-extice*dzetai_05int(i)*l1)
  enddo

  if (Ti1(1) > Meltpnt(0.d0) + T_phase_threshold) then
    dhwhigh = (Ti1(1) - Meltpnt(0.d0) - T_phase_threshold) * &
    & ci_m_roi*l1*ddzi(1)*0.5d0/(row0_m_Lwi)
    dhif = -dhwhigh*row0_d_roi
  else
    dhwhigh = 0
    dhif = 0
  endif
  dhw = snmelt*dt + dhwhigh
  dhi = dhif
  dhi0 = dhi
  ice=1;snow=1;water=0;deepice=0

  call SNOW_COND_HEAT_COEF()
  call SOIL_COND_HEAT_COEF()
  call S_DIFF(dt)
  call T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat, extice, fetch, dt)
  call SOILFORLAKE( &
  & ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat,extice,fetch,dt,a,b,c,d,Temp)
  call SNOWTEMP(ix,iy,nx,ny,year,month,day,hour,snowmass, &
  & snowmass_init,a,b,c,d,Temp,phi,extwat,extice,fetch,dt)
  call METHANE &
  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil,TgrAnn, &
  & ddz, Tw2, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
  & fplant, febul, fdiff, ftot, fdiff_lake_surf, &
  & plant_sum,bull_sum,oxid_sum,rprod_sum, &
  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC,h_talik,tot_ice_meth_bubbles,add_to_winter)
  qmethane(1:M+1) = qwater(1:M+1,2)
  qmethane(M+2:M+ns) = qsoil(2:ns)
!  call METHANE2 & ! two-meth
!  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
!  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC, & ! two-meth
!  & h_talik,tot_ice_meth_bubbles2) ! two-meth
!  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth

  l2 = l1 + dhi 
  h2 = h1 + dhw
  ls2 = 0.

else

!CASE 3.2: ICE WITHOUT SNOW !

!Radiation fluxes
  Erad = shortwave*(1-albedoofice)
  do i = 1, Mice
    SRi(i) = Erad*dexp(-extice*dzetai_05int(i)*l1)
  enddo

  if (precip>0.and. tempair<0) then
    flag_snow = 1
    hs1 = 0.02
    dhip = 0
  else
    dhip = precip*row0_d_roi*dt
  end if
  if (Ti1(1) > Meltpnt(0.d0) + T_phase_threshold) then
    dhw = (Ti1(1) - Meltpnt(0.d0) - T_phase_threshold) * &
    & ci_m_roi*l1*ddzi(1)*0.5d0/(row0_m_Lwi)
    dhihigh = - dhw*row0_d_roi
    dhis = 0
  else
    dhis = - Elatent/(roi*Liv)*dt
    dhw = 0.
    dhihigh = 0
  endif
  dhi = dhis + dhihigh + dhip
  dhi0 = dhi
  ice=1;snow=0;water=0;deepice=0 

  call SOIL_COND_HEAT_COEF()
  call S_DIFF(dt)
  call T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat, extice, fetch, dt)
  call SOILFORLAKE( &
  & ix,iy,nx,ny,year,month,day,hour,phi, &
  & extwat,extice,fetch,dt,a,b,c,d,Temp)
  call METHANE &
  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil,TgrAnn, &
  & ddz, Tw2, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
  & fplant, febul, fdiff, ftot, fdiff_lake_surf, &
  & plant_sum,bull_sum,oxid_sum,rprod_sum, &
  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC,h_talik,tot_ice_meth_bubbles,add_to_winter)
  qmethane(1:M+1) = qwater(1:M+1,2)
  qmethane(M+2:M+ns) = qsoil(2:ns)
!  call METHANE2 & ! two-meth
!  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
!  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!  & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC, & ! two-meth
!  & h_talik,tot_ice_meth_bubbles2) ! two-meth
!  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
    
  l2 = l1 + dhi
  h2 = h1 + dhw
  ls2 = 0.
 
endif

ENDIF if3 

! CASE 4: SNOW AND SOIL WITHOUT ICE AND WATER
    
if4: IF (layer_case==4) THEN
  print*, 'The non-operational case: &
  &there is no water and ice layers: STOP'
  STOP
  if (flag_snow==1) then
    ice=0;snow=1;water=0;deepice=0
    call SNOW_COND_HEAT_COEF()
    call SOIL_COND_HEAT_COEF()
    call S_DIFF(dt)
    call T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
    & extwat, extice, fetch, dt)
    call SOILFORLAKE( &
    & ix,iy,nx,ny,year,month,day,hour,phi, &
    & extwat,extice,fetch,dt,a,b,c,d,Temp)
    call SNOWTEMP(ix,iy,nx,ny,year,month,day,hour,snowmass, &
    & snowmass_init,a,b,c,d,Temp,phi,extwat,extice,fetch,dt)
    call METHANE &
    & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil,TgrAnn, &
    & ddz, Tw2, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
    & fplant, febul, fdiff, ftot, fdiff_lake_surf, &
    & plant_sum,bull_sum,oxid_sum,rprod_sum, &
    & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC,h_talik,tot_ice_meth_bubbles,add_to_winter)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns)
!    call METHANE2 & ! two-meth
!    & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!    & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
!    & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!    & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!    & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC, & ! two-meth
!    & h_talik,tot_ice_meth_bubbles2) ! two-meth
!    qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!    qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth    
  else
    ice=0;snow=0;water=0;deepice=0
    call SOIL_COND_HEAT_COEF()
    call S_DIFF(dt)
    call T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
    & extwat, extice, fetch, dt)
    call SOILFORLAKE( &
    & ix,iy,nx,ny,year,month,day,hour,phi, &
    & extwat,extice,fetch,dt,a,b,c,d,Temp)
    call METHANE &
    & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil,TgrAnn, &
    & ddz, Tw2, lammeth, qwater, h1, l1, ls1, dhw, dhw0, &
    & fplant, febul, fdiff, ftot, fdiff_lake_surf, &
    & plant_sum,bull_sum,oxid_sum,rprod_sum, &
    & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC,h_talik,tot_ice_meth_bubbles,add_to_winter)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns)
!    call METHANE2 & ! two-meth
!    & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
!    & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
!    & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
!    & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
!    & anox,M,ns,dt,rprod_total_oldC,rprod_total_newC, & ! two-meth
!    & h_talik,tot_ice_meth_bubbles2) ! two-meth
!    qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
!    qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
  endif 
ENDIF if4


!TRIBUTARY INFLOW to a lake

h2 = h2 + trib_inflow/year_sec*dt 

if (h2 < 0.005.and.ls2/=0.and.l2==0) then
  ls2 = ls2 + h2*row0_d_roi
  if (ls2 < 0.005) then 
    ls2 = 0; Tis2 = 0
  endif
  h2 = 0.
  Tw2 = 0.
  Sal2 = 0.
endif

if ((h2<0.005.and.h2>0).and.l2 /= 0) then
  l2 = l2 + h2*row0_d_roi
  if (l2 < 0.005) then 
    l2 = 0; Ti2 = 0
  endif
  h2 = 0.
  Tw2 = 0.
  Sal2 = 0. 
end if
    
if (h2 < 0.005.and.l2/=0.and.ls2/=0) then
  Ti2 = Ti2*l2/(l2+ls2) + Ti2*ls2/(l2+ls2)
  l2 = l2 + h2*row0_d_roi + ls2 
  ls2 = 0.
  Tis2 = 0.
  h2 = 0 
  Tw2 = 0.
  Sal2 = 0.
endif
   
if (l2 < 0.005.and.h2 /= 0) then
  h2 = h2 + l2*roi/row0
  if (h2 < 0.005) then 
    h2 = 0; Tw2 = 0; Sal2 = 0.
  endif
  l2 = 0.
  Ti2 = 0.
end if

if (h2 < 0.005.and. l2 == 0) then
  h2 = 0.
  Tw2 = 0.
  Sal2 = 0. 
end if

if (l2 < 0.005.and. h2 == 0) then
  l2 = 0.
  Ti2 = 0.
end if
   
if ((hs1 < 0.005.and.hs1>0).or.(hs1>0.and.l2==0.and.h2/=0)) then
! h2 = h2 + (totalprecips-totalevaps-totalmelts)
  h2 = h2 + snowmass/row0
  if (h2 < 0.005) then
    l2 = l2 + row0*h2/roi
    h2 = 0; Tw2 = 0; Sal2 = 0.
    if (l2 < 0.005) then 
      l2 = 0; Ti2 = 0
    endif
  endif
  hs1 = 0.
  flag_snow = 0
  flag_snow_init = 1
end if

if (hs1==0 .and. flag_snow==1) then
  flag_snow = 0
  flag_snow_init = 1
endif
 
if (ls2 < 0.005.and.h2/=0) then
  h2 = h2 + ls2*roi/row0
  if (h2 < 0.005) then 
    h2 = 0; Tw2 = 0; Sal2 = 0.
  endif
  ls2 = 0
  Tis2 = 0.
endif 

if (ls2>0.and.h2 == 0) then
  l2 = ls2
  Ti2 = Tis2
  if (l2 < 0.005) then 
    l2 = 0; Ti2 = 0
  endif
  ls2 = 0
  Tis2 = 0.
endif

h1    = h2
l1    = l2
ls1   = ls2
Tw1   = Tw2
Ti1   = Ti2
Tis1  = Tis2
Sal1  = Sal2
Sals1 = Sals2
qwater(:,1) = qwater(:,2)
!qwater2(:,1,1) = qwater2(:,2,1) ! two-meth
!qwater2(:,1,2) = qwater2(:,2,2) ! two-meth
oxyg  (:,1) = oxyg  (:,2)
carbdi(:,1) = carbdi(:,2)
u1    = u2
v1    = v2
Tskin(1) = Tskin(2)
    
! CALCULATION OF SUMMARY FLUXES !
totalpen = totalpen + dhwfsoil 
totalevap = totalevap + Elatent/(row0*Lwv)*dt
totalprecip = totalprecip + precip*dt
totalhflux = totalhflux + hflux*dt
totalerad = totalerad + erad*dt

if (l1==0) tsw = Tw1(1) + 273.15
if (l1/=0 .and. flag_snow==0) tsw = Ti1(1) + 273.15
if (flag_snow==1) tsw = T(itop) + 273.15

if (init(ix,iy)==0) then
  init(ix,iy)=1  
endif

!VALUES AT NEXT TIME STEP (time + dt) IN CURRENT POINT (ix,iy)

call UPDATE_NEXT_TIMESTEP ( &
& ix, iy, nx, ny, M, Mice, ns, ms, ml, &
& l1, h1, ls1, hs1, &
& u1, v1, &
& E1, eps1, k_turb_T_flux, &
& Tsoil1, Sals1, wi1, wl1, &
& Tw1, Sal1, lamw, &
& Tskin(1), &
& Ti1, Tis1, &
& dz, T, wl, dens, &
& qwater(1,1), qsoil, &
& oxyg(1,1), carbdi(1,1), &
& snmelt, &
& cdmw, &
& time, &
& dhwfsoil, &
& Elatent, &
& dhw, &
& dhw0, &
& dhi, &
& dhi0, &
& dls0, &
& velfrict, &
& roughness, &
& eflux0_kinem, &
& tot_ice_meth_bubbles, &

& flag_snow, flag_snow_init, &
& itop, &
& nstep )

hw1 = hw
xlew1 = xlew
cdmw1 = cdmw

!if (error_cov==1) then
! if (int(nstep*dt/(interval*60*60))>=nstep_meas) then
!  nstep_meas = nstep_meas + 1
!  if (nstep_meas<=ntim) then
!   do i = 1, M+1
!    j = (i-1)*n_modvar
!    x_mod(nstep_meas,j+1)= Tw2(i)
!    x_mod(nstep_meas,j+2)= E2(i)
!   enddo 
!  endif
! endif
!endif 
!if (error_cov==1.and.nstep_meas==ntim) then
! call HL_method
! error_cov = 0
!endif

!Diagnoctics

!call fluxes

!if (ix==10.and.iy==10) print*, 'Lake', &
!& 'T1=',Tw1(1),'T2=',Tw1(2),'h1=',h1,'S=',shortwave, &
!& 'hh=',hflux,'LE=',Elatent, &
!& 'eflux=',eflux, 'Longwave=', longwave, &
!& 'Bal=',(shortwave*(1-albedoofwater)*(sabs)+longwave- &
!&  surfrad-hflux-Elatent-eflux)/(cw_m_row0*ddz(1)*h1/2)*dt

if (accum == .false. .and. &
&   year  == year_accum_begin  .and. &
&   month == month_accum_begin .and. &
&   day   == day_accum_begin   .and. &
&   hour  >= hour_accum_begin) accum = .true.

if (accum == .true. .and. &
&   year  == year_accum_end  .and. &
&   month == month_accum_end .and. &
&   day   == day_accum_end   .and. &
&   hour  >= hour_accum_end) accum = .false.

if (accum) then
! Note: variables accumulation is implemented currently only for one-point simulations!
  call ACCUM_VAR &
  & (dt, l1, febul, fdiff, fdiff_lake_surf, &
  & rprod_total_newC, rprod_total_oldC, &
  & febultot, fdifftot, fdiff_lake_surftot, &
  & rprod_total_newC_integr, rprod_total_oldC_integr, &
  & add_to_winter)
!  call ACCUM_VAR & ! two-meth
!  & (dt, l1, febul2, fdiff2, fdiff_lake_surf2, & ! two-meth
!  & rprod_total_newC, rprod_total_oldC, & ! two-meth
!  & febultot2, fdifftot2, fdiff_lake_surftot2, & ! two-meth
!  & rprod_total_newC_integr2, rprod_total_oldC_integr2) ! two-meth
endif

if (flag_print) then ! The output in ASCII files 
  if (monthly_out ==1) call MON_OUT(ix,iy,nx,ny,year,month,day,hour)
  if (daily_out   ==1) call DAY_OUT(ix,iy,nx,ny,year,month,day,hour)
  if (hourly_out  ==1) call HOUR_OUT(ix,iy,nx,ny,year,month,day,hour)
  if (everystep   > 0) call EVERYSTEP_OUT(ix,iy,nx,ny)
  if (time_series == 1) then
    call SERIES_OUT(ix,iy,nx,ny,year,month,day,hour,tsw)
    ndec = 4 
    i = 1
    call PROFILE_OUTPUT &
    & (ix, iy, nx, ny, &
    &  year, month, day, hour, &
    &  time, dt_out, &
    &  Tw1, z_full, M+1, &
    &  zgrid_out, ngrid_out, &
    &  outpath, 'water_temp', i, ndec, .false.)
    ndec = 4
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, nx, ny, &
    &  year, month, day, hour, &
    &  time, dt_out, &
    &  Tsoil3, zsoil, ns, &
    &  zgridsoil_out, -1, & !ngridsoil_out, &
    &  outpath, 'soil_temp', i, ndec, .false.)
    z_watersoil(1:M+1) = z_full(1:M+1)    
    do j = 2, ns
      z_watersoil(M+j) = h1 + zsoil(j)
    enddo
    ndec = 10
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, nx, ny, &
    &  year, month, day, hour, &
    &  time, dt_out, &
    &  qmethane, z_watersoil, M+ns, &
    &  zgridsoil_out, -1, &
    &  outpath, 'methane_water_soil', i, ndec, .false.)
    ndec = 4
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, nx, ny, &
    &  year, month, day, hour, &
    &  time, dt_out, &
    &  oxyg(1,1), z_full, M+1, &
    &  zgrid_out, -1, &
    &  outpath, 'oxygen_water', i, ndec, .false.)
    ndec = 4
    i = i + 2    
    call PROFILE_OUTPUT &
    & (ix, iy, nx, ny, &
    &  year, month, day, hour, &
    &  time, dt_out, &
    &  wl1, zsoil, ns, &
    &  zgrid_out, -1, &
    &  outpath, 'wl_soil', i, ndec, .false.)
    ndec = 4
    i = i + 2    
    call PROFILE_OUTPUT &
    & (ix, iy, nx, ny, &
    &  year, month, day, hour, &
    &  time, dt_out, &
    &  wi1, zsoil, ns, &
    &  zgrid_out, -1, &
    &  outpath, 'wi_soil', i, ndec, .true.)    
  endif
endif
 
if (runmode == 1)  then
  if (mod(nstep,nscreen)==0) then
    write (*,'(3a25)') 'timestep', 'N of point', &
    & 'Surface temperature'
    write (*,'(2i25,f25.1)') nstep, ix, tsw
    write (*,'(4a8)') 'Year', 'Month', 'Day', 'Hour'
    i = year 
    if (year<1000) i = year+1000
    write (*,'(3i8,f8.2)') i, month, &
    & day   , hour
    write(*,*) 'Total summer ebul,  Total summer bottom diff,  Total summer surface diff'
    write(*,*) febultot(1), fdifftot(1), fdiff_lake_surftot(1)
    write(*,*) 'Total winter ebul,  Total winter bottom diff,  Total winter surface diff'
    write(*,*) febultot(2), fdifftot(2), fdiff_lake_surftot(2)
    write(*,*) 'Total summer meth-prod from "new" org,  Total summer meth-prod from "old" org'
    write(*,*) rprod_total_newC_integr(1), rprod_total_oldC_integr(1)
    write(*,*) 'Total winter meth-prod from "new" org,  Total winter meth-prod from "old" org'
    write(*,*) rprod_total_newC_integr(2), rprod_total_oldC_integr(2)    
!    write(*,*) 'Total summer ebul,  Total summer bottom diff,  Total summer surface diff' ! two-meth
!    write(*,*) 'old', febultot2(1,1), fdifftot2(1,1), fdiff_lake_surftot2(1,1) ! two-meth
!    write(*,*) 'new', febultot2(1,2), fdifftot2(1,2), fdiff_lake_surftot2(1,2) ! two-meth
!    write(*,*) 'Total winter ebul,  Total winter bottom diff,  Total winter surface diff' ! two-meth
!    write(*,*) 'old', febultot2(2,1), fdifftot2(2,1), fdiff_lake_surftot2(2,1) ! two-meth
!    write(*,*) 'new', febultot2(2,2), fdifftot2(2,2), fdiff_lake_surftot2(2,2) ! two-meth
!    write(*,*) 'Total summer meth-prod from "new" org,  Total summer meth-prod from "old" org' ! two-meth
!    write(*,*) rprod_total_newC_integr(1), rprod_total_oldC_integr(1) ! two-meth
!    write(*,*) 'Total winter meth-prod from "new" org,  Total winter meth-prod from "old" org' ! two-meth
!    write(*,*) rprod_total_newC_integr(2), rprod_total_oldC_integr(2) ! two-meth
  endif
endif
  
!FORMATS!
7   format (f7.3, 4i5,37f7.3) 
60  format (f5.2, 2e16.5, 2f8.3, e15.5, f11.5, f7.1)
! 61  format (i5, f6.2, 2e16.5, 2f8.3, e15.5, f11.5, f7.1)
61  format (f6.2, f5.2, 2e16.5, 2f8.3, e15.5, f11.5, f7.1)
62  format (f6.2, f5.2, 2f12.3, 2f8.3, f11.5, f7.1,3f12.4,2f7.2,3f7.1)
80  format (f7.3, 10f9.2)
90  format (f7.3, 12f9.2)   
100 format (a5, 2a16, 2a8, a15, a11, a7) 
!FINISHING PROGRAM! 
  
END SUBROUTINE LAKE

!-----------------------------------------------------------------------------------
!TEMPERATURE EQUATION SOLVER
!-----------------------------------------------------------------------------------

SUBROUTINE T_SOLVER(ix,iy,nx,ny,year,month,day,hour,phi, &
& extwat, extice, fetch, dt)

!T_SOLVER implements iterations to find water surface temperature 
!at the next time step    

use DRIVING_PARAMS
use ARRAYS
use INOUT_PARAMETERS, only : &
& lake_subr_unit_min, &
& lake_subr_unit_max

implicit none

!Input variables
integer(4), intent(in) :: ix, iy
integer(4), intent(in) :: nx, ny

integer(4), intent(in) :: year, month, day

real(8), intent(in) :: hour
real(8), intent(in) :: phi
real(8), intent(in) :: extwat, extice
real(8), intent(in) :: fetch
real(8), intent(in) :: dt

!Local variables
real(8), allocatable :: Tsurf_prev(:,:)
integer(4), allocatable :: surftyp_prev(:,:)
integer(4), allocatable :: surftyp(:,:)

real(8) :: dt_scan
real(8) :: Tsurf1, Tsurf2, Tsurf3
real(8) :: Bal1, Bal2, Bal3
real(8) :: total_Bal = 0.
real(8) :: resid_Bal_max = 0.01d0

integer(4) :: iter
integer(4) :: maxiter
integer(4) :: n_unit = lake_subr_unit_min

logical :: firstcall = .true.

!External functions
real(8), external :: HEATBALANCE

SAVE

if (firstcall) then
  call CHECK_UNIT(lake_subr_unit_min,lake_subr_unit_max,n_unit)
  open (n_unit,file=path(1:len_trim(path))// &
  & 'results/debug/iter_T.dat')
  maxiter = 10
  dt_scan = 5.

  allocate (Tsurf_prev(nx,ny))
  allocate (surftyp_prev(nx,ny))
  allocate (surftyp(nx,ny))
endif

if (snow == 1) then
  surftyp(ix,iy) = snow_indic
elseif (ice == 1) then
  surftyp(ix,iy) = ice_indic
elseif (water == 1) then
  surftyp(ix,iy) = water_indic
else
  surftyp(ix,iy) = soil_indic
endif
call SURF_CHAR(surftyp(ix,iy),year,month,day,hour,phi,fetch)
 
if (nstep > 1.and.surftyp(ix,iy) == surftyp_prev(ix,iy)) then
  Tsurf2 = Tsurf_prev(ix,iy) - 10.
elseif (surftyp(ix,iy) == water_indic) then
  Tsurf2 = -10.
else
  Tsurf2 = -90.
endif

!print*, 'In T_SOLVER: ', Tsurf2, surftyp(ix,iy)

!SCANNING INTERVAL OF SURFACE TEMPERATURE [-90 C, ...]

Bal1=1.; Bal2=1.
do while (Bal1*Bal2>0)
  Tsurf1 = Tsurf2
  Tsurf2 = Tsurf2 + dt_scan
  call T_diff(1,Tsurf1,dt)
  Bal1 = HEATBALANCE(Tsurf1,surftyp(ix,iy),extwat,extice,fetch,dt) 
  call T_diff(1,Tsurf2,dt)
  Bal2 = HEATBALANCE(Tsurf2,surftyp(ix,iy),extwat,extice,fetch,dt) 
  if (Tsurf2>80.) then
    print*, 'Severe: iterations for the surface temperature &
    &do not converge: at the point', ix, iy, 'STOP', 'nstep = ', nstep, lamw
!    print*, 'Temperature limit in scan process is exceeded: STOP'
    STOP
  endif
enddo

iter=0
Bal3=1.

!CHORDE METHOD TO FIND SURFACE TEMPERATURE
if (Bal1 /= 0.d0.and.Bal2 /= 0.d0) then
  cy1:do while (dabs(Bal3) > resid_Bal_max)
    iter = iter + 1
    call T_diff(1,Tsurf1,dt)
    Bal1 = HEATBALANCE(Tsurf1,surftyp(ix,iy),extwat,extice,fetch,dt) 
    call T_diff(1,Tsurf2,dt)
    Bal2 = HEATBALANCE(Tsurf2,surftyp(ix,iy),extwat,extice,fetch,dt) 
    Tsurf3 = (Tsurf1*Bal2 - Tsurf2*Bal1)/(Bal2 - Bal1)
    call T_diff(1,Tsurf3,dt)
    Bal3 = HEATBALANCE(Tsurf3,surftyp(ix,iy),extwat,extice,fetch,dt)
    if (Bal1*Bal3 < 0.) then
      Tsurf2 = Tsurf3
    elseif (Bal2*Bal3 < 0.) then
      Tsurf1 = Tsurf3
    endif
    if (iter > maxiter) then
      write (n_unit,*) Bal3
      exit cy1
    endif
  enddo cy1
elseif (Bal1 == 0.d0) then
  Tsurf3 = Tsurf1
elseif (Bal2 == 0.d0) then
  Tsurf3 = Tsurf2
endif

call T_diff(0,Tsurf3,dt)

Tsurf_prev(ix,iy) = Tsurf3
surftyp_prev(ix,iy) = surftyp(ix,iy)

total_Bal = total_Bal + Bal3

if (firstcall) firstcall = .false.
END SUBROUTINE T_SOLVER


SUBROUTINE T_DIFF(surf,Tsurf,dt)

!T_DIFF calculates temperature profile in water, soil, ice and snow,
!with known temparature at the surface 

use NUMERIC_PARAMS
use PHYS_CONSTANTS
use DRIVING_PARAMS
use ARRAYS
use PHYS_FUNC, only: &
& MELTPNT 
implicit none

integer(4), parameter :: vector_length = 350

!Variables
real(8) dt,Tsurf,Hflow
real(8), dimension(1:ms) :: Tsn,lams,q,cs
real(8), dimension(1:vector_length) :: a,b,c,d,Temp
real(8) AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens
real(8) dz
integer(4) i,j,itop,surf
   

common /snow_char/ Tsn,cs
common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
& ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML),dens(ms)
common /SOILDAT/ dz(ms),itop
common /watericesnowarr/ lams,q
!data num /1/

SAVE

Hflow=0.
!Meltpnt - Melting point temperature, C degrees 
!Meltpnt=0.
!ddz=1./float(M)VS,06.2007

if (surf==1) then
  if (snow==1) then
    c(itop)=1.
    b(itop)=0.
    d(itop)=Tsurf
    if (itop<=ms-2) then
      call DIFF_COEF(a,b,c,d,itop+1,ms-1,itop+1,ms-1,snow_indic,dt) 
    endif
    if (ice==1) then
!----------------SNOW-ICE INTERFACE
      a(ms)=-(lams(ms-1)+lams(ms))/(2.*dz(ms-1))
      b(ms)=-lami/(ddzi(1)*l1)+ci_m_roi*dhi0/(2.*dt)
      c(ms)=a(ms)+b(ms)-(cs(ms)*dens(ms)*dz(ms)/(2.*dt) + &
      & ci_m_roi*ddzi(1)*l1/(2.*dt))
      d(ms)=-Ti1(1)*(cs(ms)*dens(ms)*dz(ms)/(2.*dt) + &
      & ci_m_roi*ddzi(1)*l1/(2.*dt))-SR_botsnow + SRi(1)
!-----------------------------------
      call DIFF_COEF(a,b,c,d,2,Mice,ms+1,ms+Mice-1,ice_indic,dt) 
      if (water==1) then
        a(ms+Mice)=0.
        c(ms+Mice)=1.
        d(ms+Mice)=Meltpnt(Sal2(1))
        call PROGONKA (a,b,c,d,Temp,itop,ms+Mice)
        do i=itop, ms
          Tsn(i)=Temp(i)
        enddo
        do i=ms, ms+Mice
          Ti2(i-ms+1)=Temp(i)
        enddo 
      else
!----------------ICE-SOIL INTERFACE-------------------------
        a(ms+Mice)=-lami/(ddzi(Mice)*l1)+ci_m_roi*(dhi-dhi0)/(2.*dt)
        b(ms+Mice)=-(lamsoil(1)+lamsoil(2))/(2.*dzs(1))
        c(ms+Mice)=a(ms+Mice)+b(ms+Mice) - &
        & (csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
        & ci_m_roi*ddzi(Mice)*l1/(2.*dt))
        d(ms+Mice)=-Ti1(Mice+1)*(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
        & ci_m_roi*ddzi(Mice)*l1/(2.*dt))-SRi(Mice)
!----------------------------------------------------------- 
        call DIFF_COEF(a,b,c,d,2,ns-1,ms+Mice+1,ms+Mice+ns-2,soil_indic,dt)
        c(ms+Mice+ns-1)=1.
        a(ms+Mice+ns-1)=1.
        d(ms+Mice+ns-1)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
        call PROGONKA (a,b,c,d,Temp,itop,ms+Mice+ns-1)
        do i=itop, ms
          Tsn(i)=Temp(i)
        enddo
        do i=ms, ms+Mice
          Ti2(i-ms+1)=Temp(i)
        enddo 
        do i=ms+Mice, ms+Mice+ns-1
          Tsoil2(i-ms-Mice+1)=Temp(i)
        enddo
      endif
    else
!-------------------SNOW-SOIL INTERFACE-------------------------
      a(ms)=-(lams(ms-1)+lams(ms))/(2.*dz(ms-1))
      b(ms)=-(lamsoil(1)+lamsoil(2))/(2*dzs(1))
      c(ms)=a(ms)+b(ms)-(cs(ms)*dens(ms)*dz(ms)/(2.*dt) + &
      & csoil(1)*rosoil(1)*dzs(1)/(2.*dt))
      d(ms)=-Tsoil1(1)*(cs(ms)*dens(ms)*dz(ms)/(2.*dt) + &
      & csoil(1)*rosoil(1)*dzs(1)/(2.*dt))-SR_botsnow
!---------------------------------------------------------------
      call DIFF_COEF(a,b,c,d,2,ns-1,ms+1,ms+ns-2,soil_indic,dt) 
      c(ms+ns-1)=1.
      a(ms+ns-1)=1.
      d(ms+ns-1)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
      call PROGONKA (a,b,c,d,Temp,itop,ms+ns-1)
      do i=itop, ms
        Tsn(i)=Temp(i)
      enddo
      do i=ms, ms+ns-1
        Tsoil2(i-ms+1)=Temp(i)
      enddo
    endif
  elseif (ice==1) then
    b(1)=0.
    c(1)=1.
    d(1)=Tsurf
    call DIFF_COEF(a,b,c,d,2,Mice,2,Mice,ice_indic,dt)  
    if (water==1) then
      a(Mice+1)=0.
      c(Mice+1)=1.
      d(Mice+1)=Meltpnt(Sal2(1))
      call PROGONKA (a,b,c,d,Temp,1,Mice+1)
      do i=1,Mice+1
        Ti2(i)=Temp(i)
      enddo
    else
!----------------ICE-SOIL INTERFACE-------------------------
      a(Mice+1)=-lami/(ddzi(Mice)*l1)+ci_m_roi*(dhi-dhi0)/(2.*dt)
      b(Mice+1)=-(lamsoil(1)+lamsoil(2))/(2*dzs(1))
      c(Mice+1)=a(Mice+1) + b(Mice+1) - &
      & (csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
      & ci_m_roi*ddzi(Mice)*l1/(2.*dt))
      d(Mice+1)=-Ti1(Mice+1)*(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
      & ci_m_roi*ddzi(Mice)*l1/(2.*dt))-SRi(Mice)
!-----------------------------------------------------------
      call DIFF_COEF(a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt) 
      c(Mice+ns)=1.
      a(Mice+ns)=1.
      d(Mice+ns)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
      call PROGONKA (a,b,c,d,Temp,1,Mice+ns)
      do i=1,Mice+1
        Ti2(i)=Temp(i)
      enddo
      do i=Mice+1,Mice+ns
        Tsoil2(i-Mice)=Temp(i)
      enddo
    endif
  elseif (water==1) then  
    b(1)=0.
    c(1)=1.
    d(1)=Tsurf
    call DIFF_COEF(a,b,c,d,2,M,2,M,water_indic,dt) 
    if (deepice==1) then
      a(M+1)=0.
      c(M+1)=1.
      d(M+1)=Meltpnt(Sal2(M+1))
      call PROGONKA(a,b,c,d,Temp,1,M+1) 
      do i=1,M+1
        Tw2(i)=Temp(i)
      enddo
    else
!---------------------WATER-SOIL INTERFACE------------------
      a(M+1)=-lamw(M)/(ddz(M)*h1)+cw_m_row0*(dhw-dhw0)/(2.*dt) + &
      & 0.5d0*cw_m_row0*PEMF(M)
      b(M+1)=-(lamsoil(1)+lamsoil(2))/(2*dzs(1))
      c(M+1)=a(M+1)+b(M+1)-(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
      & cw_m_row0*ddz(M)*h1/(2.*dt))-cw_m_row0*PEMF(M)
      d(M+1)=-Tw1(M+1)*(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
      & cw_m_row0*ddz(M)*h1/(2.*dt))-SR(M) - &
      & cw_m_row0*PEMF(M)*pt_down_f(M)
!-----------------------------------------------------------
      call DIFF_COEF(a,b,c,d,2,ns-1,M+2,M+ns-1,soil_indic,dt)
      c(M+ns)=1.
      a(M+ns)=1.
      d(M+ns)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
      call PROGONKA (a,b,c,d,Temp,1,M+ns)
      do i=1,M+1
        Tw2(i)=Temp(i)
      enddo
      do i=M+1,M+ns
        Tsoil2(i-M)=Temp(i)
      enddo
    endif
  else
    b(1)=0.
    c(1)=1.
    d(1)=Tsurf  
    call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_indic,dt) 
    c(ns)=1.
    a(ns)=1.
    d(ns)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
    call PROGONKA (a,b,c,d,Temp,1,ns)
    do i=1,ns
      Tsoil2(i)=Temp(i)
    enddo
  endif
else   
  if (water==0) then
    RETURN
  else
    if (ice==1) then
      b(1)=0.
      c(1)=1.
      d(1)=Meltpnt(Sal2(1))
      call DIFF_COEF(a,b,c,d,2,M,2,M,water_indic,dt)  
      if (deepice==1) then
        a(M+1)=0.
        c(M+1)=1.
        d(M+1)=Meltpnt(Sal2(M+1))
        call PROGONKA (a,b,c,d,Temp,1,M+1)
        do i=1,M+1
          Tw2(i)=Temp(i)
        enddo
        b(1)=0.
        c(1)=1.
        d(1)=Meltpnt(Sal2(M+1))
        call DIFF_COEF(a,b,c,d,2,Mice,2,Mice,deepice_indic,dt)
!----------------DEEPICE-SOIL INTERFACE-------------------------
        a(Mice+1)=-lami/(ddzi(Mice)*ls1)+ci_m_roi*(dls-dls0)/(2.*dt)
        b(Mice+1)=-(lamsoil(1)+lamsoil(2))/(2*dzs(1))
        c(Mice+1)=a(Mice+1)+b(Mice+1) - &
        & (csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
        & ci_m_roi*ddzi(Mice)*ls1/(2.*dt))
        d(Mice+1)=-Tis1(Mice+1)*(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
        & ci_m_roi*ddzi(Mice)*ls1/(2.*dt))-SRdi(Mice)
!-----------------------------------------------------------
        call DIFF_COEF(a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt)  
        c(Mice+ns)=1.
        a(Mice+ns)=1.
        d(Mice+ns)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
        call PROGONKA (a,b,c,d,Temp,1,Mice+ns)
        do i=1,Mice+1
          Tis2(i)=Temp(i)
        enddo
        do i=Mice+1,Mice+ns
          Tsoil2(i-Mice)=Temp(i)
        enddo
      else
!---------------------WATER-SOIL INTERFACE------------------
        a(M+1)=-lamw(M)/(ddz(M)*h1)+cw_m_row0*(dhw-dhw0)/(2.*dt) + &
        & 0.5d0*cw_m_row0*PEMF(M)
        b(M+1)=-(lamsoil(1)+lamsoil(2))/(2*dzs(1))
        c(M+1)=a(M+1)+b(M+1)-(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
        & cw_m_row0*ddz(M)*h1/(2.*dt))-cw_m_row0*PEMF(M)
        d(M+1)=-Tw1(M+1)*(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
        & cw_m_row0*ddz(M)*h1/(2.*dt))-SR(M) - &
        & cw_m_row0*PEMF(M)*pt_down_f(M)
!-----------------------------------------------------------
        call DIFF_COEF(a,b,c,d,2,ns-1,M+2,M+ns-1,soil_indic,dt)
        c(M+ns)=1.
        a(M+ns)=1.
        d(M+ns)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
        call PROGONKA (a,b,c,d,Temp,1,M+ns)
        do i=1,M+1
          Tw2(i)=Temp(i)
        enddo
        do i=M+1,M+ns
          Tsoil2(i-M)=Temp(i)
        enddo
      endif
    elseif (deepice==1) then    
      b(1)=0.
      c(1)=1.
      d(1)=Meltpnt(Sal2(M+1))
      call DIFF_COEF(a,b,c,d,2,Mice,2,Mice,deepice_indic,dt)
!----------------DEEPICE-SOIL INTERFACE-------------------------
      a(Mice+1)=-lami/(ddzi(Mice)*ls1)+ci_m_roi*(dls-dls0)/(2.*dt)
      b(Mice+1)=-(lamsoil(1)+lamsoil(2))/(2*dzs(1))
      c(Mice+1)=a(Mice+1)+b(Mice+1) - &
      & (csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
      & ci_m_roi*ddzi(Mice)*ls1/(2.*dt))
      d(Mice+1)=-Tis1(Mice+1)*(csoil(1)*rosoil(1)*dzs(1)/(2.*dt) + &
      & ci_m_roi*ddzi(Mice)*ls1/(2.*dt))-SRdi(Mice)
!-----------------------------------------------------------
      call DIFF_COEF(a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt)  
      c(Mice+ns)=1.
      a(Mice+ns)=1.
      d(Mice+ns)=Hflow*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
      call PROGONKA (a,b,c,d,Temp,1,Mice+ns)
      do i=1,Mice+1
        Tis2(i)=Temp(i)
      enddo
      do i=Mice+1,Mice+ns
        Tsoil2(i-Mice)=Temp(i)
      enddo
    endif
  endif
endif

END SUBROUTINE T_DIFF
    



!---------------------------------------------------------------------------------


SUBROUTINE DIFF_COEF(a,b,c,d,n0,n1,m0,m1,substr,dt)

!-------------------DEFINES COEFFICIENTS FOR SOLVING A THERMAL DIFFUSIVITY-----------------
!-------------------EQUATION BY FACTORIXATION METHOD----------------------------------------

use NUMERIC_PARAMS
use PHYS_CONSTANTS  
use DRIVING_PARAMS
use ARRAYS
implicit none

integer(4), parameter :: vector_length = 350

!Input variables
real(8), intent(in) :: dt

integer(4), intent(in) :: n0
integer(4), intent(in) :: n1
integer(4), intent(in) :: m0
integer(4), intent(in) :: m1
integer(4), intent(in) :: substr

real(8), intent(out), dimension(1:vector_length) :: a,b,c,d

!Common variables
integer(4) :: itop
real(8), dimension(1:ms) :: Tsn, cs
real(8), dimension(1:ms) :: lams, q
real(8) :: dz
real(8) :: AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens

common /snow_char/ Tsn,cs
common /watericesnowarr/ lams,q
common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
& ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML),dens(ms)
common /SOILDAT/ dz(ms),itop

!External functions
real(8), external :: DZETA
real(8), external :: DZETAI

!Local variables
real(8) :: lam1
real(8) :: lam2
real(8) :: dzmean

integer(4) :: i
integer(4) :: j

SAVE

if (m1-m0/=n1-n0) then
  print*, 'Error: diff_coef'
  STOP
endif

SELECT CASE (substr)
!CASE 1: SOIL
  case (soil_indic)
    do i = m0, m1
      j = i - m0 + n0
      lam2 = 0.5d0*(lamsoil(j)+lamsoil(j+1))
      lam1 = 0.5d0*(lamsoil(j-1)+lamsoil(j))
      dzmean = 0.5d0*(dzs(j-1)+dzs(j))
      a(i)=-lam1/dzs(j-1)
      b(i)=-lam2/dzs(j)
      c(i)=-(lam1/dzs(j-1)+lam2/dzs(j)) - &
      & dt_inv*(csoil(j)*rosoil(j)*dzmean)
      d(i)=-Tsoil1(j)*dt_inv*(csoil(j)*rosoil(j)*dzmean)
    enddo 
!CASE 2: ICE
  case (ice_indic)
    do i = m0, m1
      j = i - m0 + n0
      c(i) = -(lami/ddzi(j-1)+lami/ddzi(j))/l1 - &
      & dt_inv*ddzi05(j-1)*l1*ci_m_roi
      a(i) = -lami/(ddzi(j-1)*l1) + &
      & ci_m_roi*dhi*dzetai_int(j)*dt_inv05 - &
      & ci_m_roi*dhi0*dt_inv05
      b(i) = -lami/(ddzi(j)*l1) - &
      & ci_m_roi*dhi*dzetai_int(j)*dt_inv05 + &
      & ci_m_roi*dhi0*dt_inv05  
      d(i) = -Ti1(j)*dt_inv*ddzi05(j-1)*l1*ci_m_roi + &
      & SRi(j) - SRi(j-1) 
    enddo
!CASE   3: SNOW  
  case (snow_indic)
    do i = m0, m1
      j = i - m0 + n0
      lam2 = 0.5d0*(lams(j)+lams(j+1))
      lam1 = 0.5d0*(lams(j-1)+lams(j))
      dzmean = 0.5d0*(dz(j-1)+dz(j))
      a(i) = -lam1/dz(j-1)
      b(i) = -lam2/dz(j)
      c(i) = -(lam1/dz(j-1)+lam2/dz(j)) - &
      & dt_inv*(cs(j)*dens(j)*dzmean)
      d(i) = -T(j)*dt_inv*(cs(j)*dens(j)*dzmean)
    enddo 
!CASE 4: WATER
  case (water_indic)
    do i = m0, m1
      j = i - m0 + n0
      c(i) = - (area_half(j)*lamw(j)/ddz(j) + &
      & area_half(j-1)*lamw(j-1)/ddz(j-1))/(h1*area_int(j)) - &
      & dt_inv*ddz05(j-1)*h1*cw_m_row0 + &
      & cw_m_row0*0.5d0*(PEMF(j)-PEMF(j-1))
      a(i) = - area_half(j-1)*lamw(j-1)/(ddz(j-1)*h1*area_int(j)) + &
      & cw_m_row0*dhw*dzeta_int(j)*dt_inv05 - &
      & cw_m_row0*dhw0*dt_inv05 + cw_m_row0*0.5d0*PEMF(j-1)
      b(i) = - area_half(j)*lamw(j)/(ddz(j)*h1*area_int(j)) - &
      & cw_m_row0*dhw*dzeta_int(j)*dt_inv05 + &
      & cw_m_row0*dhw0*dt_inv05 - cw_m_row0*0.5d0*PEMF(j)
      d(i) = - Tw1(j)*dt_inv*ddz05(j-1)*h1*cw_m_row0 + &
      & SR(j) - SR(j-1) - cw_m_row0*PEMF(j-1)*pt_down_f(j-1) + &
      & cw_m_row0*PEMF(j)*pt_down_f(j)
    enddo 
!CASE 5: DEEP ICE
  case (deepice_indic)
    do i = m0, m1
      j = i - m0 + n0
      c(i) = -(lami/ddzi(j)+lami/ddzi(j-1))/ls1 - &
      & dt_inv*ddzi05(j-1)*ls1*ci_m_roi
      a(i) = -lami/(ddzi(j-1)*ls1)+ci_m_roi*dls*dzetai_int(j)* &
      & dt_inv05-ci_m_roi*dls0*dt_inv05
      b(i) = -lami/(ddzi(j)*ls1)-ci_m_roi*dls*dzetai_int(j)* &
      & dt_inv05+ci_m_roi*dls0*dt_inv05
      d(i) = -Tis1(j)*dt_inv*ddzi05(j-1)*ls1*ci_m_roi + &
      & SRdi(j) - SRdi(j-1)
    enddo
!CASE 6: SALINITY IN WATER
 case (water_salinity_indic)
   do i = m0, m1
     j = i - m0 + n0
     c(i)=-(lamsal(j)/ddz(j)+lamsal(j-1)/ddz(j-1))/h1 - &
     & dt_inv*ddz05(j-1)*h1  
     a(i)=-lamsal(j-1)/(ddz(j-1)*h1)+dhw*dzeta_int(j)*dt_inv05 - &
     & dhw0*dt_inv05
     b(i)=-lamsal(j)/(ddz(j)*h1)-dhw*dzeta_int(j)*dt_inv05 + &
     & dhw0*dt_inv05 
     d(i)=-Sal1(j)*dt_inv*ddz05(j-1)*h1   
   enddo 
!CASE 7: SALINITY IN THE SOIL  
 case (soil_salinity_indic)
   do i = m0, m1
     j = i - m0 + n0
     a(i) = -wsoil(j-1)*dt/(dzs(j-1)+dzs(j)) 
     b(i) = wsoil(j)*dt/(dzs(j-1)+dzs(j))   
     c(i) = -1.d0-wsoil(j)*dt/(dzs(j-1)+dzs(j)) + &
     & wsoil(j-1)*dt/(dzs(j-1)+dzs(j))   
     d(i) = -Sals1(j)  
   enddo 
!CASE 8: METHANE IN WATER     
 case (water_methane_indic)
   do i = m0, m1
     j = i - m0 + n0
     c(i)=-(lammeth(j)/ddz(j)+lammeth(j-1)/ddz(j-1))/h1 - &
     & dt_inv*ddz05(j-1)*h1  
     a(i)=-lammeth(j-1)/(ddz(j-1)*h1)+dhw*dzeta_int(j)*dt_inv05 - &
     & dhw0*dt_inv05
     b(i)=-lammeth(j)/(ddz(j)*h1)-dhw*dzeta_int(j)*dt_inv05 + &
     & dhw0*dt_inv05
     d(i)=-qwater(j,1)*dt_inv*ddz05(j-1)*h1   
   enddo
!CASE 9: OXYGEN IN WATER
 case (water_oxygen_indic)
   do i = m0, m1
     j = i - m0 + n0
     c(i) = - (lamoxyg(j)/ddz(j)+lamoxyg(j-1)/ddz(j-1))/h1 - &
     & dt_inv*ddz05(j-1)*h1  
     a(i) = - lamoxyg(j-1)/(ddz(j-1)*h1)+dhw*dzeta_int(j)*dt_inv05 - &
     & dhw0*dt_inv05
     b(i) = - lamoxyg(j)/(ddz(j)*h1)-dhw*dzeta_int(j)*dt_inv05 + &
     & dhw0*dt_inv05 
     d(i) = - oxyg(j,1)*dt_inv*ddz05(j-1)*h1   
   enddo
!CASE 10: CARBON DIOXIDE IN WATER
 case (water_carbdi_indic)
   do i = m0, m1
     j = i - m0 + n0
     c(i) = - (lamcarbdi(j)/ddz(j)+lamcarbdi(j-1)/ddz(j-1))/h1 - &
     & dt_inv*ddz05(j-1)*h1  
     a(i) = - lamcarbdi(j-1)/(ddz(j-1)*h1)+dhw*dzeta_int(j)*dt_inv05 - &
     & dhw0*dt_inv05
     b(i) = - lamcarbdi(j)/(ddz(j)*h1)-dhw*dzeta_int(j)*dt_inv05 + &
     & dhw0*dt_inv05 
     d(i) = - carbdi(j,1)*dt_inv*ddz05(j-1)*h1   
   enddo          
END SELECT
  
END SUBROUTINE DIFF_COEF


SUBROUTINE SURF_CHAR(surftyp,year,month,day,hour,phi,fetch)
use DRIVING_PARAMS
use PHYS_CONSTANTS
use ARRAYS
use ATMOS, only: &
& WIND, &
& VELFRICT_PREV
use SFCFLX, only: &
& SFCFLX_ROUGHNESS    
use PHYS_FUNC, only: &
& SINH0, &
& WATER_ALBEDO, &
& CHARNOCK_Z0 

!----------------DEFINES SURFACE CHARACTERISTICS:-------------------
!----------------roughness,emissivity,albedo,relative humidity,-----
!----------------coefficients in Magnus formula---------------------


implicit none

integer(4), intent(in) :: surftyp
integer(4), intent(in) :: year
integer(4), intent(in) :: month
integer(4), intent(in) :: day
real(8)   , intent(in) :: hour
real(8)   , intent(in) :: phi
real(8)   , intent(in) :: fetch

real(8) roughness,emissivity,albedo,aM,bM,relhums
real(8) x1,x2,x3,x4,x5

common /surface/ roughness,emissivity,albedo,aM,bM,relhums
 
SAVE 
    
select case (surftyp)
  case (1)
    albedo = albedoofsoil
    roughness = 0.05
    aM = aMagw 
    bM = bMagw
    emissivity = emissivityofsoil
    relhums = 0.5
  case (2) 
    albedo = albedoofice
    emissivity = emissivityofice
    aM = aMagi
    bM = bMagi 
    roughness = 0.00001
    relhums = 0.7
  case (3)
    albedo = albedoofsnow
    emissivity = emissivityofsnow
    aM = aMagi
    bM = bMagi 
    roughness = 0.001
    relhums = 0.7
  case (4)
    velfrict_prev = dmax1(velfrict_prev, 1.d-2)
!   Extended Charnock formula
    roughness = CHARNOCK_Z0(velfrict_prev)
!   call SfcFlx_roughness (fetch, wind, velfrict_prev, 0.d0, &
!   & x1, x2,roughness, x4, x5)
!   roughness = 1.d-3
    aM = aMagw
    bM = bMagw
    if (varalb==0) then
      albedo = albedoofwater
    elseif (varalb==1) then
      albedo = WATER_ALBEDO( SINH0(year,month,day,hour,phi) )
      !albedo = (dirdif()*0.05/(sinh0()+0.15)+0.05)/ &
      !(1+dirdif())
    endif
    emissivity = emissivityofwater 
    relhums = 1. !0.9
end select  

END SUBROUTINE SURF_CHAR



FUNCTION DZETA(nl)
use ARRAYS, only : ddz
use DRIVING_PARAMS, only : M

!The function DZETA
!returns the dzeta coordinate [0..1]
!(non-dimensional coordinate in water layer)
!of the level nl [1..M+1].
!nl might be fractional number, however the function is opimized for
!integer or half-integer nl.

implicit none

real(8) :: DZETA

!Input variable
real(4), intent(in) :: nl

!local variables
real(8), allocatable, save :: dzeta_int(:)
real(8), allocatable, save :: dzeta_05int(:)

integer(4) :: nl_int
integer(4) :: i !Loop index

logical :: firstcall
data firstcall /.true./

if (firstcall) then
  allocate (dzeta_05int(1:M))
  dzeta_05int(1) = 0.5d0*ddz(1)
  do i = 2, M
    dzeta_05int(i) = dzeta_05int(i-1) + 0.5d0*(ddz(i-1) + ddz(i))
  enddo
  allocate (dzeta_int(1:M+1))
  dzeta_int(1) = 0.d0
  do i = 2, M+1
    dzeta_int(i) = dzeta_int(i-1) + ddz(i-1)
  enddo
endif

if (mod(nl,1._4) == 0._4) then
  DZETA = dzeta_int(int(nl))
elseif (mod(2._4*nl,1._4) == 0._4) then
  DZETA = dzeta_05int(int(nl))
else
  nl_int = int(nl)
  if (nl_int > 1) then
    DZETA = sum(ddz(1:nl_int-1))
  else
    DZETA = 0.
  endif
  if (nl<M+1) DZETA = DZETA + (nl-real(nl_int))*ddz(nl_int)
endif

if (firstcall) firstcall = .false.
END FUNCTION DZETA


FUNCTION DZETAI(nl)
use ARRAYS, only : ddzi
use DRIVING_PARAMS, only : Mice

!The function DZETAI
!returns the dzetai coordinate [0..1]
!(non-dimensional vertical coordinate in ice layers)
!of the level nl [1..M+1].
!nl might be fractional number, however the function is optimized
!for integer or half-integer nl.

implicit none

real(8) :: DZETAI

real(4), intent(in) :: nl

real(8), allocatable, save :: dzetai_int(:)
real(8), allocatable, save :: dzetai_05int(:)

integer(4) :: nl_int
integer(4) :: i !Loop index

logical :: firstcall
data firstcall /.true./

if (firstcall) then
  allocate (dzetai_05int(1:Mice))
  dzetai_05int(1) = 0.5d0*ddzi(1)
  do i = 2, Mice
    dzetai_05int(i) = dzetai_05int(i-1) + 0.5d0*(ddzi(i-1) + ddzi(i))
  enddo
  allocate (dzetai_int(1:Mice+1))
  dzetai_int(1) = 0.d0
  do i = 2, Mice + 1
    dzetai_int(i) = dzetai_int(i-1) + ddzi(i-1)
  enddo
endif

if (mod(nl,1._4) == 0._4) then
  DZETAI = dzetai_int(int(nl))
elseif (mod(2._4*nl,1._4) == 0._4) then
  DZETAI = dzetai_05int(int(nl))
else
  nl_int = int(nl)
  if (nl_int>1) then
    DZETAI = sum(ddzi(1:nl_int-1))
  else
    DZETAI = 0.
  endif
  if (nl<Mice+1) DZETAI = DZETAI + (nl-real(nl_int))*ddzi(nl_int)
endif

if (firstcall) firstcall = .false.
END FUNCTION DZETAI


REAL(8) FUNCTION VARMEAN(var,grid)
use arrays
use driving_params

!The function varmean computes
!the spatial average of the
!variable var.
!If grid = 1, than var is averaged as 
!if it has been placed in integer points of mesh  
!If grid = 2, than var is averaged as 
!if it has been placed in half-integer points of mesh  

implicit none

integer(4) grid,i
real(8), intent(in):: var(1:M+1)

varmean = 0.
if(grid==1) then
  varmean = varmean + var(1)  *0.5*ddz(1)
  varmean = varmean + var(M+1)*0.5*ddz(M)
  do i=2,M
    varmean = varmean + var(i)*0.5*(ddz(i-1)+ddz(i))
  enddo
elseif (grid==2) then
  do i=1,M
    varmean = varmean + var(i)*ddz(i)
  enddo
else
  print*, 'Illegal identifier GRID in the function varmean: STOP'
  STOP
endif

RETURN
END FUNCTION VARMEAN
